Modeling gravitational radiation from coalescing binary black holes 
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With the goal of bringing theory, particularly numerical relativity, to bear on an astrophysical 
problem of critical interest to gravitational wave observers we introduce a model for coalescence 
radiation from binary black hole systems. We build our model using the Lazarus approach, a 
technique that bridges far and close limit approaches with full numerical relativity to solve Einstein 
equations applied in the truly nonlinear dynamical regime. We specifically study the post-orbital 
radiation from a system of equal-mass non-spinning black holes, deriving waveforms which indicate 
strongly circularly polarized radiation of roughly 3% of the system's total energy and 12% of its 
total angular momentum in just a few cycles. Supporting this result we first establish the reliability 
of the late-time part of our model, including the numerical relativity and close-limit components, 
with a thorough study of waveforms from a sequence of black hole configurations varying from 
previously treated head-on collisions to representative target for "ISCO" data corresponding to the 
end of the inspiral period. We then complete our model with a simple treatment for the early 
part of the spacetime based on a standard family of initial data for binary black holes in circular 
orbit. A detailed analysis shows strong robustness in the results as the initial separation of the 
black holes is increased from 5.0 to 7.8M supporting our waveforms as a suitable basic description 
of the astrophysical radiation from this system. Finally, a simple fitting of the plunge waveforms 
is introduced as a first attempt to facilitate the task of analyzing data from gravitational wave 
detectors. 

PACS numbers: 04.25.Dm, 04.25. Nx, 04.30.Db, 04.70. Bw 



I. INTRODUCTION 

Black holes may be the most interesting and extraor- 
dinary objects predicted by Einstein's general relativity. 
Over the last forty years, a flowering of observational 
techniques and theoretical consideration has produced 
a change in perspective, bringing black holes from sus- 
pected mathematical fictions to the presumptive pow- 
erhouses behind the most energetic astronomical phe- 
nomena under study today. From quasars and Gamma 
Ray Bursts (GRB's) to ordinary galactic nuclei and the 
final state of sufficiently large stars, black holes seem 
now quite common astrophysical objects. The charac- 
teristic of black holes which has brought them to such 
prominence in astrophysics is that when they are not 
isolated they are emphatically not black but produce the 
most luminous astrophysical phenomena known. One can 
estimate that a coalescing binary system of two black 
holes releases energy with a peak luminosity of about 
lO^'^c^/G, 10^'^ times the power output of the Sun, dur- 
ing the very brief period of its merger. 

These events have not yet been observed because, lack- 
ing matter, the binary black hole (BBH) system releases 
its energy purely in the form of gravitational waves. 
To the developing field of gravitational wave astronomy 
these are expected to be generally the brightest sources 
in the sky. For ground-based gravitational observato- 
ries, such as LIGO [§ and GEO600, followed shortly by 
VIRGO M and several resonant bar detectors, all sen- 



sitive to "higher frequencies" in the 10^ - lO'^ Hz band, 
the relevant coalescence events have masses on the or- 
der of large stars lO'^ - 10^ Mq. Such systems may be 
generated through the evolution of stellar binaries, or 
as recent simulations suggest through many-body inter- 
actions in globular clusters^, Q]. Space based detectors, 
such as the proposed NASA-ESA LISA mission [|[ 0] will 
be sensitive to a lower frequency band roughly 10~ - 10*^ 
Hz and the relevant systems are super-massive black hole 
binaries involving black holes in galactic cores. 

For both classes of observation, the importance of the- 
oretical models for the expected BBH coalescence is ex- 
traordinary ^ . For the current generation of ground- 
based interferometers beginning to take scientific data 
in the next few months, the anticipated signal-to- noise 
ratios are sufficiently small that theoretical information 
from BBH coalescence models may be essential to distin- 
guishing genuine detections from noise events. For LISA, 
on the other hand, the signal-to-noise ratios for super- 
massive BBH events should be quite high, up to 10'* 
In this case the BBH events may be so strong and nu- 
merous that the ability to make observations of other 
weaker sources depends on filtering out the BBH signals 
according to accurate model waveforms. Because the new 
ground-based interferometers are already beginning op- 
erations and because model information may be critical 
even in the developmental stages of the LISA mission, 
there is a pressing need for producing a fair model for 
the BBH coalescence immediately. As general relativis- 
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tic modeling techniques continue to evolve and improve, 
there is naturally a trade-ofF between the level of confi- 
dence we can have in a model and how long we wait to 
produce it, but the clear and present need for some kind 
of theoretical result, provides a strong argument against 
waiting. Our goal here is to introduce a practical first 
model which may be progressively improved as new tech- 
niques are developed. Producing coalescence waveforms 
from the first usable models opens up a vital channel of 
communication between observational and theoretical ef- 
forts that is now crucially required. Of course we expect 
this theoretical information to support the observational 
effort, providing observers with specific examples of what 
sort of signal they can expect. Importantly, though less 
widely recognized in the theoretical community, a provi- 
sional model allows observational practicalities to inform 
the theoretical effort. Not all theoretical concerns are 
equally consequential in practice, and enhanced interac- 
tion with observers through model waveforms will help to 
identify which details of the model need the most work. 

Clearly the foundation for a theoretical model of the 
BBH system will be Einstein's general theory of rela- 
tivity. Intensive efforts to develop numerical codes able 
to solve the Einstein field equations using state-of-the- 
art supercomputers have been underway for more than a 
decade now . These efforts have resulted in successful 
studies of several phenomena but, the treatment of bi- 
nary black hole systems has proven very difficult. Thus 
far, the full numerical treatment of these systems in three 
dimensions (3D), based on the traditional 34-1 decompo- 
sition of space and time, has been severely limited on one 
hand, because of huge computer memory and processing 
requirements, and, on the other, by instabilities related 
to the formulation of the evolution equations. These dif- 
ficulties combine to cause codes to fail to be accurate 
before any useful gravitational wave information can be 
directly extracted even for sample problems such as, dis- 
torted black holes, and the so-called non-axisymmetric 
'grazing collisions', where the black holes must start out 
yet too close to allow a clear astrophysical interpretation 
Ipr] , |l2| (and even for single distorted black holes @). 

In an earlier letter we presented the first theoret- 
ical estimates for the gravitational radiation waveforms 
and energy to be expected from the plunge of orbiting 
non-spinning binary black holes starting from an esti- 
mate of the innermost stable circular orbit (ISCO) as 
provided by Cook/Baumgarte jl^, |l6), using effective po- 
tential method applied to Bowen-York initial data. Such 
a binary system has also been considered by Buonanno 
and Damour , who produce a radiation waveform es- 
timate without solving the full system of Einstein's equa- 
tions in the (expected) nonlinear regime, but using ex- 
trapolations of the 2.5 Post-Newtonian order. In this 
paper we apply the Lazarus approach p8[ |, to develop a 
functional astrophysical model for binary black hole co- 
alescence. 

A complete implementation of our approach requires 
at least five components. There are three basic treat- 



ment regimes: (i) a Far-Limit (FL) treatment covering 
the earliest region of spacetime when the black holes are 
still relatively far apart, (ii) a Fully-Nonlinear (FN) ap- 
proach to the essentially nonlinear plunge/merger region 
of the spacetime and (iii) a Close-Limit (CL) approach 
to the dynamics of the newly formed final black hole and 
the propagation of radiation. In addition to these are 
two interfaces, FL-FN and FN-CL, which require a spe- 
cific means of propagating the relevant physical informa- 
tion from the end of one regime to the beginning of the 
next. Section II provides an overview of our model and 
a brief review of the FN, FN-CL, CL techniques which 
are described in detail and tested in principle in Ref. 
]l8[ . In Section [II C we apply theses techniques to per- 



form a systematic study of the waveforms generated by 
a set of sample problems, simulations of equal mass bi- 
nary black holes with no intrinsic spin starting from near 
the ISCO. This provides a proving ground for testing our 
machinery for taking initial data at the beginning of the 
FN regime and producing a final waveform result. The 
study connects to previously well-studied problems and 
builds up to the practical demands of our astrophysical 
model, with Spinning and unequal mass black holes can 
also be treated within our method and will be published 
in a different paper. 

To complete the model we need to provide FL and 
FL-FN treatments. As our approach to this regime has 
not been treated extensively in Ref. jl^ we provide here 
a more detailed exposition. In Sec. Ill we describe our 
treatment and discuss how this approach compares with 
other candidate treatments which might also seem appro- 
priate concluding that based on kinematics ours is the 
best approach to providing modeled astrophysical data 
to feed into into the FN part of the calculation available 
which is so far sufficiently developed for practical appli- 
cations. But a specific advantage of having a complete 
model is that we can move beyond these (often dubious) 
kinematical arguments and explicitly study the robust- 
ness of the full dynamical model. We begin this analysis 
in Sec. with successful robustness test of the FL and 
FL-FN part of the model. We also extend our resul ts to 
meet Post-Newtonian (PN) calculations in Sec. VP . We 
evolve black holes from the Post-Newtonian determined 
parameters for the ISCO. Section VI presents summarizes 



the astrophysically relevant results of our treatment for 
the coalescence waveforms, and provides a simple analyt- 
ical approximation of our results as a practical represen- 
tation for the benefit gravitational wave observers and 
data analysts. 



II. THE MODEL 

Under the Lazarus approach our model is specified by 
providing a specific implementation for five components 
(FL, FL-FN, FN, FN-CL, CL). A sketch of this com- 
bined approach can be outlined by the following steps: 
(1) First provide description of the early dynamics of the 
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system with a FL approach, such as the PN method, 
which is appropriate for slowly moving, well-separated 
black holes, or some alternative quasi-stationary (QS) 
method. (2) Extract critical information about the late- 
time configuration of this system in terms of Cauchy 
data, Qij and Kij, and translate this information to a 
corresponding solution of the gravitational initial-value 
problem. (3) Apply a full 3D numerical simulation of 
Einstein's equations to generate a numerical spacetime 
covering the non-linear interaction region. The evolution 
should proceed for long enough so that the subsequent 
evolution of the region exterior to the final single rem- 
nant black hole can be well approximated by perturba- 
tive dynamics. (4) At this point we choose a "late-time" 
slice from the numerically generated spacetime, extract 
■04 = CajS-ysn^fh^n^fh^ and dt4'4, to quantify the de- 
viation of the numerical spacetime from a Kerr geome- 
try. Then (5) evolve via the Teukolsky equation, which 
governs the dynamics of Kerr perturbations in the time- 
domain p9| , long enough to drive all significant radiation 
into the radiation zone where it can be interpreted. 

For convenience in discussion we will refer to steps 
(3-5) through which radiation waveforms are derived 
from numerical simulation initial data as the "late-stage" 
part of the model. Likewise the "early-stage" part of 
the model encompasses steps (1-2) in which numerical 
data are produced from some model of the early part of 
the spacetime. A suitable implementation of the late- 
stage part of the model has been presented in Ref. [ p^ . 
During its development, the late-stage implementation 
has passed many tests, including the robustness of ip4 
and dtipi against nonlinear contamination of the Cauchy 
data In a first letter our method proved to be 
capable of determining radiation waveforms for a model 
problem, the head-on collision of two black holes, with 
accuracy comparable to the best published 2D numer- 
ical results, allowing at the same time a more direct 
physical understanding of the collisions and indicating 
clearly when non-linear dynamics are important as the 
final black hole is formed. In | |l^ we report on internal 
consistency checks after including net angular momen- 
tum in the system, e.g. quadratic convergence to vanish- 
ing gravitational radiation, which our method passes, for 
the evolution of Kerr initial data. This is a non-trivial 
test of our procedure since the computed spurious radia- 
tion energy and waveforms will give us a direct measure 
of the 'internal' error with which we can determine such 
quantities. Typically, the levels of spurious radiation we 
have found for rotation parameter a/M = 0.8 is around 
10-^Af. 



Below, in Sec. Ill C we supplement these tests with a 
developmental study of the performance of our late-stage 
implementation on a series of test problems which build 
from the previous test problems to a fiducial case of prac- 
tical interest, evolution from "ISCO" data, as described 
below. 

In this work we will add to the model a simple treat- 
ment for the the early-stage part. We essentially model 



the early part of the spacetime as a stack of spacelike 
slices, each slice chosen from a family of solutions of the 
Einstein initial value problem representing non-spinning, 
orbiting black holes. The solutions on the slices are fur- 
ther selected by kinematical (effective potential) argu- 
ments to correspond to circular orbits. This early-stage 
part of our model is clearly less developed than what 
we have implemented for the late-stage, and it is fair to 
note at the start that there are many reasons to pursue 
improvements for this part of the model. The key advan- 
tage of our present approach is that it can be applied even 
without further development to provide a complete pro- 
visional model which produces coalescence waveforms. It 
is precisely in the context of such a completed model that 
we have a concrete basis for e valuating any treatment of 
steps (1-2). In Section IV A we provide a detailed de- 



scription of our treatment and introduce some relevant 
concerns related to this part of the model. In particular 
we consider the plausibility of our model in comparison 
with PN calculations. 

One of the most attractive characteristics of the 
Lazarus approach is that it provides a built-in facility 
for cross-checking the three primary treatments underly- 
ing the model. By shifting the location of the interface 
regions in the model spacetime, we can effectively ex- 
change the treatment chosen for that region of spacetime. 
This creates an inherent facility for validating the three 
evolutionary model components. We can cross-check the 
Close Limit perturbation theory treatment with Numer- 
ical Relativity results for the same region of spacetime, 
systematically incorporating a common technique for ver- 
ifying both Close Limit and numerical results |]l3|, |2^, ^ . 
In |]l8| we described the basis for such comparisons via 
the FN-CL interface looking for phase locking and en- 
ergy plateaux in the waveforms. With the addition of 
the early-stage treatment we can now begin to practically 
compare the FL and FN components of the model. Sec- 
tion ^ is dedicated to evaluating the FL model by com- 
parison with numerical simulation. Although, the model 
for the FL is clearly more suspect, comparison of the 
results are in much better agreement than expected, val- 
idating our simple treatment as a fair provisional model. 



III. PRACTICAL TESTING OF THE 
LATE-STAGE MODEL 



Our idea for the late-stage part of the model is to use 
the Lazarus tools for deriving radiation waveforms, en- 
ergy and momentum radiated, plunge durations, etc, re- 
sulting from an evolution of binary black holes that have 
started close 'enough' to the location of the ISCO. We 
approach this case developmentally, beginning with an 
axisymmetric system of ISCO-separation black holes re- 
leased from rest. This system has been studied exten- 
sively using 2D numerical simulations and with a pre- 
vious simplified version of our model's late-stage treat- 
ment pl|. From this head-on collision case we gradu- 
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ally add orbital angular momentum to the system to ap- 
proach the target ISCO problem. We call this sequence of 
initial data sets the 'P-sequence' through which we vary 
the transverse momentum of the individual black hole 
from zero, in the head-on case to our target P — Pisco- 
The waveforms we present here are the culmination of 
a number of numerical simulations which were needed 
to establish both the effectiveness of our methodology 
and the plausibility of these results as a fair astrophysi- 
cal model. The first two subsections lay the groundwork 
with a description of the particular family of initial data 
we are studying, and an overview of our methodology. 



P=0.33 M 



=0.45 M 



L = 4.9M 



m=0.45M 



J = 0.77 M 



P=0.33M 



FIG. 1: Effective potential - ISCO data 



A. Initial Data 

The Cauchy data required to begin a 3-f 1 simulation 
of vacuum general relativity can be given in the form of a 
metric tensor g^j on an initial spacelike and an extrinsic 
curvature tensor Kij which contains information about 
how the initial hypersurface is embedded in the space- 
time. The gravitational field equations include four con- 
straints, limiting the choices for gij and Kij. To present 
data for a binary black hole system generally requires 
(1) some astrophysical ansatz for plausible field configu- 
rations around (and constituting) the two black holes (2) 
solving the constraints. 

Many conceivable approaches have been proposed and 
to varying degrees developed, for the case of black hole 
initial data. In this study we apply a variant of the 
best known approach, known as Bowen-York data [ p4[ . 
This approach provides an ansatz for describing black 
hole systems through a set of restrictions on the met- 
ric gij and extrinsic Kij. The assumptions are: con- 
formal flatness gij = ^'^6ij, a maximal slice Kl = 0, 
and a purely longitudinal form for the extrinsic curva- 
ture K'^ $-iO(V'V"^ + V^Y' - Ig'^VkV"). Instead 
of the traditional "inversion-symmetric" boundary con- 
ditions at each of the black holes, we apply a simplified 
variant ('puncture' data [^) which allows black hole 
data to be represented numerically without excising the 
interior regions of the black holes. An advantage of this 
class of data is that is has been studied extensively. 

A significant feature, expected to exist in the binary in- 
teraction, is the ISCO indicating the closest possible con- 
figuration of black holes at which the conservative part 
of the dynamics allow stable orbits. Although the field 
theory nature of the 'Full Numerical' approach does not 
lend itself naturally to a particle like description of the 
dynamics implicit in describing the ISCO, it is possible to 
construct particle-like parameterizations of appropriately 
chosen families of solutions to the initial value problem. 
In these terms, the kinematical relationships, among the 
solutions to the initial value problem provide a descrip- 
tion of the system similar to the conservative particle 
dynamics, including an estimate for the ISCO. 

Different approaches have been developed to identify 
quasi-circular orbital configurations and to determine the 



location and frequency of the ISCO. We will use here the 
results based on the effective potential method of [^sj as 
derived in JT6[ for the puncture construction of black hole 
initial data. This provides a space of solutions parame- 
terized by the momenta, spins and positions of each black 
hole. From these solutions one can apply the Arnowitt- 
Deser-Misner (ADM) formalism to determine the total 
energy Mad m for the slice of spacetime. If a bare mass 
is inferred for each individual black hole from its appar- 
ent horizon area, then the difference EEff can be used 
to identify points in the solution space corresponding to 
quasi-circular orbits, and ISCO. 

For puncture data the ISCO is characterized by the 
parameters 



L/M = 4.9, P/M = 0.335, 
J/M^ = 0.77, m = 0.45M. 



(3.1) 



where m is the mass of each of the single black holes, 
M is the total ADM mass of the binary system, L is 
the proper distance between the apparent horizons, P is 
the magnitude of the linear momenta (equal but opposite 
and perpendicular to the line connecting the holes) and 
J is the total angular momentum. This configuration is 
represented in Fig. |l| 

Families of solutions to the initial value problem have 
been produced from other sets of assumptions on the 
fields as well. An approach inspired from the Kerr-Schild 
form of a Kerr black hole [Em, provides a more natu- 
ral model at least when the black holes are far apart, 
but has not yet been studied as a kinematical model for 
quasi-circular orbits and ISCO. Another treatment along 
the lines of the maximal/conformally flat formalism with 
no assumption about the longitudinality of the extrin- 
sic curvature but with additional restrictions guarantee- 
ing that these conditions are initially consistent with the 
evolution equations, is very promising but has not yet 
been applied to irrotational black hole systems as consid- 
ered here 28, 29j. Finally, an interesting first step in 
providing astrophysically realistic initial based FN mo- 
tivated choices for the freely-specifiable portions of the 
data jsot will be subject of further investigation within 
our approach. 
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B. methods 

A detailed discussion of our techniques has been pre- 
sented in a previous paper [T^ so we only briefly review 
the techniques here. We first perform a numerical simu- 
lations to cover a relatively brief period of fundamentally 
non-linear two-body interactions to the point where the 
dynamics of the final black hole dominate and can be 
treated linearly. We then integrate the standard ADM 
form of Einstein's equations using the thrice iterated 
Crank-Nicholson approach with maximal slicing and van- 
ishing shift, making extensive use of the publicly available 
Cactus computational toolkit and Briigmann's multi- 
grid elliptic solver BAM for initial data and maximal 
slicing choices. In designing these simulations we have 
focused on the particular needs of our problem defined 
within the Lazarus approach, rather than on alternative 
techniques for generic applications. Some of the priorities 
in this case are: optimal accuracy for brief simulations 
in the entire region outside the final horizon, avoidance 
of boundary noise, and consistency with gauge assump- 
tions applied in interfacing to the "close-limit" phase of 
the treatment. The accuracy requirement tends to fa- 
vor the standard ADM over other approaches. Although 
the ADM system tends to be more susceptible to insta- 
bilities, it also tends to be more accurate before these 
instabilities suddenly arise. 

In the cases considered here, these runs would suffer 
catastrophic instabilities after 15 — 20M evolution, but 
the more critical question for our problem is that the 
runs be sufficiently accurate in the convergent regime 
prior to the instability. In some cases, other systems 
seem to allow for significant improvements in the treat- 
ment of boundaries. But at present even these improved 
boundary treatments remain problematic, allowing phys- 
ical wave content into the computational domain. For 
brief, accurate simulations it is important, and practical 
to move the boundary far away, causally separating it 
from the relevant wave-generation region. To realize this 
without undue computational expense we cast the initial 
data in specialized fisheye coordinates which bring a dis- 
tant outer boundary to a closer coordinate location with- 
out sacrificing physical resolution in the strong field re- 
gions ijl^ . Our gauge assumptions favor maximal slicing 
in the numerical simulations, and for simplicity without 
an obvious alternative, we use a vanishing shift. Other 
sophisticated techniques currently under development in 
the field, such as black hole excision, are not so relevant 
to our narrowly defined numerical simulation needs. We 
effectively excise the interior of the final event horizon 
when we make the transition to (purely exterior) close 
limit perturbation theory. 

In analyzing our results we make extensive use of sev- 
eral special techniques. A key tool for interpreting our 
numerically simulated spacetimes is the complex-valued 
'iS-invariant ' . 

S = 27J^/I^, (3.2) 



where I and are the two complex curvature invariants 
X and JJ, which are essentially the square and cube of 
the self-dual part, Ca/i^s = C'a/j^a -I- (i/2)eafc„mC™^, of 
the Weyl tensor: 

I ^ Co^p^sC'^''''^ and J = Cc.f3.,sC''l,C''-"'^ . (3.3) 

In vacuum general relativity this is the unique (up to 
rescaling) spacetime scalar quantity which can be con- 
structed from the components of the curvature tensor 
which is not dominated by its fall-off behavior in an 
asymptotically fiat spacetime. It is thus a natural tool 
for invariant interpretation of spacetime dynamics. A key 
characteristic of the 5-invariant which we take advantage 
of here is that is has the exact value, 5 = 1 in the Kerr 
geometry and tends to provide a useful measure of devi- 
ations from the Kerr background. As a rule of thumb, 
we adopt the criterion that black hole perturbation the- 
ory is likely to be effective if the real part of S differs 
from unity by less than a factor of two (when reaching 
the maximum) outside of our estimated location for the 
horizon. In these binary black hole spacetimes, we find 
the maximal values of S on the orbital (z) axis. 

Our most important means of evaluating the calcu- 
lations is the inherent cross-checking between the two 
evolutionary spacetime models (close limit and numeri- 
cal simulation). By changing the location of the transi- 
tion between these approaches in our model spacetime, 
which is achieved by varying the evolution time T of the 
numerical simulations, we can effectively exchange one 
treatment for the other. If the resulting radiation is ap- 
proximately unchanged, then the two independent treat- 
ments are agreeing with each other and we can be confi- 
dent in the results. In practice we compare the results by 
separately examining the magnitude (total radiation en- 
ergy E) and phasing of the waveforms. Thus we look at 
the energy as a function of the transition time T, and ex- 
pect to find a plateau in E once we have evolved through 
the most significant regions of non-linear dynamics and 
before the inaccuracies in the numerical simulation have 
grown dramatically. 



C. P-Sequence results 

In order to better understand the physics of the 
plunge we have designed a set of sequences approach- 
ing the ISCO by changing one of its physical parame- 
ters. Many different sequences are possible. The 'P- 
sequence' for which we keep the separation constant at 
L = Ljsco = 4.9M, but vary the linear momentum, 
P/ Pisco = 0, 1/3, 2/3, 5/6, lis particularly interesting 
because it connects the well-studied case of head-on col- 
lisions from rest to our target, ISCO data. Fig. ^ illus- 
trates the initial configurations studied with increasing 
linear momentum from the resting holes to the ISCO 
values. 
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m=0.45 M 



J = 0.77A^ P=0.33M 
J = 0.5M^ P=0.22M 
J = 0.25 AT P=0.11M 



P=0.11 M 
P=0.22 M 
P=0.33 M 



m=0.45 M 



O 
O 



-2 -■ 1 2 3 



L = 4.9M 



FIG. 2: Sequence P 



(Brill-Lindquist initial data) 





The head-on colUsion of black holes from rest has been 
studied using our techniques in Ref. pl| ] and compared 
with 2D simulations. Here we revisit this configuration, 
now not with the holes lying along the 0— axis, but along 
the J/— axis, and having switched from Misner to Brill- 
Lindquist initial data, to serve as a reference point for 
the ISCO initial data. The change in orientation implies 
that the radiation is now predominated by is m = ±2 
component, each contributing a factor 3/8 to the total 
energy; while m = waves have similar time dependence 
but smaller, contributing to 1/4 of the total radiated en- 
ergy. Note that the reflection symmetry of the systems 
considered here implies that the m = -1-2 and m — —2 
modes are directly related and contribute the same to 
the total energy radiated. Without loss of information 
we present only the m = +2 results. 

Fig. 1^ shows four snapshots of the apparent horizon 
surfaces on the orbital plane js^. They show that the 
two black holes start well detached. Clearly, the 'grid 
stretching' effect due to the vanishing shift we used dur- 
ing evolution makes them appear to grow in the coordi- 
nate space already after 8M of full numerical evolution. 
These plots are particularly useful to extract qualitative 
information about the system. Soon after a common ap- 
parent horizon covers the system which tends to have an 
increasingly close to spherical shape at later times. This 
tends to set an upper limit on the time at which linear 
theory can take over as has been discussed in Ref. . It 
is expected that a common event horizon should have ap- 
peared several M's of time earlier. The key physical fea- 
ture which actually makes the close limit approximation 
effective is that the black holes share a common poten- 
tial barrier which appears even earlier than the common 
event horizon. 

In Fig. ^ we show the time evolution (actually com- 
puted every \M) of the 'iS-invariant'. The 'iS-invariant' 



FIG. 3: Apparent horizons for Pjd ~ 




FIG. 4: S-invariant measuring deviations from Kerr 



clearly shows that space time is asymptotically Kerr 
{S = 1) at large distances. It also shows that the geome- 
try at the initial slice has significant distortions near the 
region where the common horizon will form and hence we 
do not expect perturbation theory to be a good approx- 
imation at very early times. In fact we need at least AM 
of evolution to bring the distortions outside the horizon 
down the level of a factor 1.5 in the iS-invariant. 

Fig. 1^ shows the computed total radiated energy from 
the system using linearized theory after T of full nonlin- 
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m=0, equator, r=20M, P=0 




FIG. 5: Total radiated energy plots 



ear numerical evolution. We see that if we insist on using 
perturbation theory to compute the subsequent evolu- 
tion. Even at very early times, we obtain that the radi- 
ated energy varies notably for different transition times. 
Only when we really reach the linear regime there is a cer- 
tain 'energy plateau'. This again happens for T « 4Af 
for both relevant modes m = 2, 0. Note that the curves 
settle down to values closely bearing the 1.5 factor among 
mode contributions expected for axisymmetric I — 2 ra- 
diation at this orientation. 

Fig. |6| displays the superposition of waveforms for the 
two relevant modes from 5M to 9M of nonlinear evo- 
lution. These waveforms are unequivocally determined 
once the transition time has been chosen so there are 
no free parameters to be adjusted in this comparison. 
Throughout this section we show the m = and m — 2 
waveforms for observers respectively on the equator and 
the z-axis, where the amplitude of each component is 
maximal. Specially useful, as an independent test of lin- 
earization, is the locking of the phases observed after 5M 
of evolution. At later times we recognize the quasinormal 
frequency of a the final Schwarzschild black hole with a 



period 



qnm 



= 16.8M. 



We observe that all the four criteria used here to de- 
termine the time this system begins to behave linearly 
outside the horizon coincide in setting it around 4 — 5M 
of nonlinear evolution. Though we have not displayed it 
explicitly, we also find good agreement with our earlier 
Misner data results. 



2. P = Pisco/S 

Now we begin to add angular momentum to the sys- 
tem. In this case we give the holes an initial momentum 
perpendicular to the separation vector of one third of the 
ISCO momentum magnitude. 

We show in Fig. |^ four snapshots of the apparent hori- 
zon on the orbital plane; a picture quite similar to the 
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FIG. 6: Waveforms for Pm = 



head-on case, except, as expected, for a slight asymmetry 
along the orbital motion and the appearance of a com- 
mon apparent horizon after 9M of full nonlinear evolu- 
tion, approximately IM later than in the head-on case. 
The asymmetry is a result of the orbital-like motion with 
respect to our vanishing shift coordinate system. The 
coordinate singularity ("puncture") which must remain 
inside the apparent horizon is fixed in place in the grid 
coordinates so that the apparent horizon cannot "orbit" 
in these coordinates, but is, rather, stretched around in 
the (/)-direction. The same characteristic is evident for 
all our runs. This twisting of the spacetime is a generic 
consequence of using a vanishing shift for these simula- 
tion and illustrates the importance of using a co-rotating 
shift with some (/j— component. 

Fig. 1^ shows the deviations from a Kerr background 
along the z— axis (the direction of maximal distortion) as 
measured by the 5— invariant. The picture shows clearly 
that at least 4M of evolution are necessary to reach a 
perturbative regime. 

In Fig. H we observe that the radiated energy reaches 
a plateau after transition times T ~ 4M consistently for 
both modes m = ±2, 0. Even for this weak radiation the 
approximate constancy of the radiated energy holds up 
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FIG. 7: Apparent horizons for Pm ~ Pisco /3 
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FIG. 8: S-invariant measuring deviations from Kerr 
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to T ~ lOM. 

Fig. |l^ displays the impressive agreement among wave- 
forms for different extraction times and different modes. 
All this reached without any adjustable parameter. The 
frequency of the waveforms quite closely resemble those 
of the least damped quasinormal modes for a Kerr black 
hole with rotation parameter a/AI — 0.26, which for 
£ = 2, m = 2 has a period pSl] r„„m ~ 15.2Af . 



FIG. 10: Waveforms for Pid = Pisco /'i 

3. P = 2Pisco/3 

This case has an initial transverse linear momentum of 
two thirds the ISCO magnitude. 

Fig. |l^ displays clear effects of orbital motion deform- 
ing the apparent horizon surfaces on the equatorial plane. 



9 




-3-2-10123 -3-2-10125 



0.00000 M 10.0000 M 




-3-2-10123 -3-2-10123 



1 1.0000 M 14.0000 M 
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FIG. 12: S-invariant measuring deviations from Kerr 



A common apparent horizon appears at T ^ \\M show- 
ing that the orbital component delays the merger of the 
holes into a single one as expected. 

Fig. |lj shows that, consistently, the 5— invariant esti- 
mate for linearization is also delayed in comparison with 
the head-on or near head-on cases. At least 8M of evo- 
lution is needed to settle down the perturbations to an 
small portion of the background. 

Likewise Fig. O shows how the plateau in the energy 
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FIG. 13: Total radiated energy plots 
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FIG. 14: Waveforms for Pm = 2Piscoli 



is reached at later times after an oscillation around the 
final value. It takes 9Af — lOM of evolution to obtain a 
stable plateau. 

Fig. |l^ displays the leading mode, m = +2, along the 
z— axis with good agreement for the times where we es- 
timate linearization takes place. Excellent phase locking 
is evident for the first three cycles with good agreement 
for the rest of the relevant signal. The frequency of the 
last part of the waveform agrees with the least damped 
quasinormal mode of a Kerr hole with rotation parameter 
a/M = 0.51 for the mode m = 2 which has a period 

of Tqnm ~ 13M . 

For the three sets of simulations just reported we have 
used grids of 384^ x 192 size (Making use of the symmetry 
along the z— axis of the problem) with a inhomogeneous 
distribution of the points called 'Fish-Eye' Each run 
required up to 70Gb of RAM memory on a Hitachi SR- 
8000 at the LRZ in Garching, Germany and took about 
10 hours running on 16 nodes (128 processors). 
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FIG. 15: Apparent horizons for Pid ~ 5Pisco /6 
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4. P = ^Piscol^ 

To closely approach the ISCO parameters we included 
the case of P = ^Piscol^- 

Fig. |l5| shows how strongly the apparent horizon looks 
deformed in the grid coordinates responding to the trans- 
verse momentum. The four snapshots shows the forma- 
tion of a common apparent horizon for T^fj = 14 — 15M 
of nonhnear evolution. Consistently, Fig. |l6| shows the 
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FIG. 17: Total radiated energy plots 



oscillations of the 5— invariant around the Kerr value, 
1, and indicates that the linearization regime is reached 
after about 9Af of full numerical evolution. Note the 
stronger deformation of the initial data labeled as T = 
compared to the previous elements of the P-sequence. 

A coherent picture comes from Fig. |l^ which suggests 
that energy plateau arises after 9M of nonlinear evolu- 
tion, just before numerical error begins to noticeably add 
to the energy.. The energy radiated oscillates around 
the final value and is evidently dominated by the modes 
m = ±2, with very little contribution from the mode 
TO = and negligible from the m — ±4 ones. We note 
that adding orbital momentum in moving along this P- 
sequence dramatically enchances the m — 2 component 
of the radiation with relatively little effect on the to = 
mode. Thus while the energy content of the two compo- 
nents was roughly comparable in the head-on case, as 
we approach ISCO the to = 2 component is growing 
strongly as would be expected for a non-axisymmetric 
orbital system. A linear analysis indicates that the 
leading ^ = 2, to = 2 contributes to the radiated energy 
with a quadratic dependence in the momentum. This is 
precisely what we observe for small P (See Fig. |2l] ). 

Fig. |l^ show excellent agreement of waveforms for dif- 
ferent transition times in the linear regime. The phase 
locking is excellent up to t 60A/, then there seems 
to be some interference from spurious radiation from the 
boundaries of the full numerical simulation that specially 
affects at transition times T > 12M. 

To compute the emitted radiation and waveforms we 
have made a first correction to the Kerr background to 
take into account the portion of the energy radiated ~ 
1.3% and the angular momentum radiated ending with 
a black hole with mass parameter Mpcerr ~ 0.987 Minu 
and ttKerr ~ O-GMpcerr- The normal frequency of the 
waveforms shown in Fig. |l^ correspond quite closely to 
the least damped quasinormal modes of this final Kerr 
hole. 

This simulation and the one we will report in the next 
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FIG. 18: Waveforms for Pm = 5Pi. 
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subsection for the ISCO havei been performed on a 512^ x 
256 grid requiring 10 hours of running time on 64 nodes 
(512 processors) of the Hitachi SR-8000 at the LRZ. The 
central resolution was M/24 and the boundaries located 
at 37M where the resolution reached ^ IM. 



D. Plunge from the ISCO data 

The ISCO configuration proves to be a borderline case 
to study with the present full numerical technologies, we 
have chosen the standard ADM formulation of Einstein 
equations in order to ensure good accuracy during the rel- 
evant nonlinear evolution. It is known that in the ADM 
formulation the instabilities that kill the evolution appear 
suddenly and do not affect the early part of the numeri- 
cal integration. This is reflected in the computation 
of the norm of the Hamiltonian constraint, H 



- Norm 



Ett2 
ij k ij k 



ijk 



(3.4) 



where Nijk is the total number of grid points labeled by 
ijk, as a first measure of the numerical errors generated 
during the unconstrained evolution. 

In Fig. |l^ we show this norm for different central res- 
olutions he- Even though we see the expected quadratic 
convergence, the exponential error growth is catastrophic 
and only slightly delayed by applying significantly more 
computational resources to increase the resoultion. Even- 
tually these runs are killed by an non-convergent instabil- 
ity, but the numerical errors are already a problem earlier 
while the code still appears to perform convergently. We 
take over with perturbation theory before the error has 
grown prohibitively. 

Fig. |2^ displays a detailed (every IM) evolution of 
the real part of the iS— invariant index showing distor- 





FIG. 20: S-invariant measuring deviations from Kerr 



tions from the Kerr background (value equal 1). Initially 
the distortions are large and the subsequent oscillations 
maintain a significant amplitude near the common hori- 
zon. After the second bounce around T = IIM the 
distortions begin to stabilize. At larger radii a signal 
with relatively large amplitude begins to leave the sys- 
tem. This will generate the first burst of radiation coming 
out from the collision. 

In Fig. ^ we gather the estimated linearization times 
for the elements of the 'P-sequence' and the ISCO and 
compare them with the time for the formation of a com- 
mon apparent horizon. They seem to carry a similar 
dependence and off set by 4 — 5M of further evolution. 
The offset is expected on the grounds of relating the lin- 
earization time to the applicability of the close limit ap- 
proximations, associated to a common potential barrier, 
which appears earlier in the evolution than a common 
event horizon and even earlier than a common apparent 
horizon. 
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FIG. 21: Linearization time and radiated energy 



We have also considered two other sequences of black 
hole configurations in the near-ISCO regime. The a- 
sequence describes configurations with constant separa- 
tion and magnitude of the hnear momentum but varies 
the angle between P and the fine joining the holes. We 
have chosen a = 90°, 60°, 30°, 0° in order to consider dif- 
ferent radial components. As an independent test we 
have also studied a sequence which connects the ISCO to 
the close limit, (where numerical simulation should not 
be needed) , by varying the initial separation of the holes 
from the ISCO separation to L/Ljsco = 1/4,1/8. For 
the sake of a compact description we are not going to 
describe those results in detail apart from reporting the 
radiated energy in Fig. |l|. Note that the 'a— sequence' 
generates more radiation than the 'P-sequence' as ex- 
pected due to its extra radial component. In this plot 
we also report on the case P = 13/l2Pjsco to check the 
effects of overshooting the ISCO data, finding consistent 
results. 

To compute the total radiated energy and angular mo- 
mentum radiated as a function of the transition time T 
displayed in Fig. ^2|, we have taken into account a change 
in the background mass Mxerr = 0.97Minitiai and the ro- 
tation parameter a Kerr = 0.7 M Kerr- We have also con- 
sidered an extrapolation of the results from two different 
resolutions in the convergent part of the full nonlinear 
simulation. In this case we find that the energy plateau 
is reached a little earlier (as expected if the background is 
closer to the actual one) and remains relatively constant 
from T = 8M through T = UM. 
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FIG. 22: Energy and angular momentum radiated 
ISCO, r'=31M, z-axis, m=+2 
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FIG. 23: ISCO Waveforms 



As remarked in Ref the c omp utation of the an- 
gular momentum radiated (Eq. 5.2) is a very sensitive 
quantity depending on correlation of waveforms. So we 
can only estimate it to be around 12% of the initial mo- 
mentum. 

The waveforms of Fig. ^ represent our best knowledge 
of the plunge so far. We have included two transition 
times in the quasi linear regime to show the differences 
and as a measure of the internal errors of the method. We 
also display the waveform obtained if no nonlinear evolu- 
tion is performed, labeled as T = OM. The clear differ- 
ences with the ones obtained in the linear regime show 
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FIG. 24: ISCO spectrum 

how important is to have some nonhnear evolution that 
settles down the system to a single rotating black hole 
plus distortions that will be radiated away (and down 
onto the final black hole). 

An interesting feature of the waveforms studied is that 
they are dominated at later times by the least damped 
quasinormal modes of the final Kerr black hole. For m = 
+2 they have a period of ~ VIM. This is refiected in the 
spectrum of the total radiated energy, including the two 
polarizations and integration over all angles, presented in 
Fig. |2|. 

This spectra, again shown for different transition times 
to have a measure of the internal error, have a rather nar- 
row peak at a frequency near that of the least damped 
quasinormal mode for a Kerr hole with ajM ~ 0.7. The 
lower frequency part of the spectrum is dominated by the 
quasinormal mode m = — 2 component with a period of 
^ 19M that is close to the radiation-orbital frequency 
(one half the orbital period of the ISCO ~ 37M). For a 
binary system of total mass M = 35M0 these frequen- 
cies correspond to A75Hz and 300Hz respectively. 



E. Critical Assessment 

We have applied our late-stage model to a sequence of 
initial data sets ranging from the head-on collision case 
to the ISCO data. Having achieved the expected result 
in the head-on case we have slowly added orbital angular 
momentum. As expected the orbital motion leads to a 
dramatic increase in the leading / = 2, m = 2 component 
of the radiation, initially depending quadratically on the 
momentum of the system. Our study proves that some 
fully non-linear treatment is essential, for a reasonable 
treatment of the BBH system in the plunge regime. On 
the other hand we find also that the extent of the space- 
time which must be treated nonlinearly in this regime is 



modest. Our estimates for the rough linearization time 
from our ISCO-candidate data set come to around lOM 
of nonlinear evolution. While the overall results at the 
ISCO are only marginally robust, within say 20% rather 
than much finer precision in the head-on case, many fea- 
tures of the waveforms are precisely determined. The 
first two cycles or so, up to the peak at t — 50M or 
so are quite precisely determined in phase. In the worst 
case, that we have underestimated our linearization time, 
and our numerical accuracy, we would not expect dra- 
matic changes in the early part of the waveform from 
an improved calculation. In many cases there are some 
evidently non-robust features in the later parts of the 
waveforms. We have traced back to several sources of 
error which influence this part of the waveform, most 
notably numerical differencing error originating near the 
punctures in our simulations and radiation errors associ- 
ated with the finite spatial extent of our numerical do- 
main. This study was a vital tool in our efforts to reduce 
finite boundary effects to the point that the overall wave- 
form seems to be accurate within 10-20%. Because the 
numerical error associated with the punctures shows up 
predominately in the higher I parts of the waveform, it is 
practical to filter these effects by extracting the strongly 
dominant / — 2 component of the radiation. This filter- 
ing preserves some 90% of the m — 2 radiation energy. 
The filtered waveforms are much more robust, and in- 
deed simpler to characterize (no- longer 2D). Therefore, 
in the subsequent analysis of our astrophysical waveform 
models we will focus on the / = 2 radiation component. 



IV. AN EARLY-STAGE MODEL 

Two important classes of treatments are applicable to 
binary black hole systems toward the end of the orbital 
epoch. In addition to 'Full Numerical' approaches which 
attempt to solve Einstein's nonlinear gravitational field 
equations are FN techniques, which make perturbative 
expansions of Einstein's equations in powers of v /c, and 
are thus most effective for describing the exterior field of 
slowly moving, separated binary systems. 

Having established a connection to previous studies al- 
lows us confidence in our calculations as solutions of Ein- 
stein's equations. But solving Einstein's equations is not 
sufficient for work which aims to treat astrophysical phe- 
nomena. For this we need a good model of the system 
to be studied, and a method for evaluating the veracity 
of the model. While we continue, separately, to develop 
an interface between the FN method and numerical sim- 
ulations which is ultimately expected to be required to 
provide the best astrophysical initial data, such data are 
not yet practical. Instead, we adopt a very basic early- 
stage model which takes no information directly from 
the FN treatment. We start our simulations with initial 
data derived from an alternative description black holes 
in quasi-circular orbit derived in Refs. [ p^ , [iot for con- 
formally flat black hole initial data, which we discussed 
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in detail in Sec. Ill A. As there are a variety of com- 
parable schemes which might be applied to provide such 
data we also need a means to assess the "astrophysical 
appropriateness" of this particular scheme. Ultimately 
the resolution of these issues will be settled by perform- 
ing a sequence of evolutions beginning with more and 
more separated initial systems approaching the regime 
where different descriptions of the initial configurations 
(Post-Newtonian, quasistationary relativistic ansatz, etc) 
merge together. 

In Sec. [VA we describe such a sequence and in 
Sec. [VB| compare our results with PN calculations appli- 
cable for well-separated systems. The comparison seems 
to justify our initial model as a plausible , a priori agree- 
ing with PN analysis to the level of accuracy presently 
available. The course for future research is to refine these 
calculations pursuing ever improved initial models begin- 
ning the simulation earlier and earlier in the astrophysical 
inspiral orbital process, thereby evaluating and improv- 
ing t he as trophysical accuracy of these early results. In 
Sec. [V A we take the first few steps along this course. 
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FIG. 25; Quasi-circular orbital data 



A. The QC sequence 

Our "early-stage" model is composed essentially of a 
sequence of spacetime slices (initial data sets). We take 
these from the same family of solutions of the initial value 
equations which we have used in Sec. Ill, "puncture- 
data" , and select data corresponding approximately to 
black holes in circular orbit. We adapt Cook's results, 
applying the effective-potential method to non-spinning 
black holes [p^, to our case to identify a sequence of 
quasicircular orbits. This selection, which we call the 
"QC" sequence, includes data ranging in initial separa- 
tion from 5M, the ISCO value, to lAM. Strictly speak- 
ing Cook's results apply to a family of data determined 
with a different interior boundary condition than applied 
for our puncture-data, but the practical differences are 
very small. Baumgarte||l^ has performed an effective 
potential analysis for the puncture data showing results 
almost identical to Cook, consistent with the noted nu- 
merical similarity of the two treatments |39 . Though 
not ideal, it is not true that these conformally flat ini- 
tial data should completely fail to approximate an as- 
trophysical configuration of black holes. For sufficiently 
well-separated binaries of non-spinning black holes the 
uncertainties associated with this method can be made 
arbitrarily small, and comparisons with post-Newtonian 
calculations have shown a reasonable correspondence in 
this limit. The black hole configurations associated 
with our QC sequence are illustrated in Fig. |2^. 

Other families of initial data are available which might 
be applied for an alternative, similarly constructed 'QC- 
sequence', but none are as yet sufficiently developed for 
this application, and there is no clear way to evaluate 
the the various a priori preferences adopted in these var- 
ious methods without a dynamical study, such as that of 



TABLE I: Quasi-Circular data 



Name 


L/M 


±X/M 


±P/M 




Am 


m/M 


QC-B |l| 


4.90 


1.158 


0.334 


0.773 


0.176 


0.450 


QC-0 |§ 


4.99 


1.169 


0.333 


0.779 


0.168 


0.453 


QC-1 


5.49 


1.364 


0.286 


0.781 


0.142 


0.463 


QC-2 


5.86 


1.516 


0.258 


0.784 


0.127 


0.470 


QC-3 


6.67 


1.849 


0.2148 


0.794 


0.1019 


0.477 


QC-4 


7.84 


2.343 


0.1743 


0.817 


0.0760 


0.483 


QC-5 


8.84 


2.772 


0.1514 


0.8397 


0.0612 


0.487 


QC-6 


9.95 


3.257 


0.1332 


0.8677 


0.0497 


0.489 


QC-7 


11.11 


3.776 


0.1190 


0.8985 


0.0408 


0.492 


QC-8 


12.17 


4.251 


0.1091 


0.9270 


0.0348 


0.494 


QC-9 


13.31 


4.77 


0.1005 


0.9584 


0.0297 


0.496 


QC-10 


14.22 


5.19 


0.0947 


0.9826 


0.0267 


0.498 



Sec.^ With our simple model we make no attempt to 
contain information from the previous dynamics of the 
spacetime {i.e radiation), but attempt a fair representa- 
tion of the system's bulk motion at late times. 

We have converted the QC data into actual physical 
parameters to be used as the initial input of the full nu- 
merical simulations. In Table Q L is the proper distance 
between the throats of the holes, X is the coordinates 
of the 'punctures' in the confornial space, P is the indi- 
vidual linear momentum of the holes, J the total angular 
momentum of the configuration, f2 the angular velocity as 
measured at infinity, and m is the bare mass of the punc- 
tures. All quantities are normalized to the total ADM 
mass of the system M. Baumgarte's determination, QC- 
B, of the ISCO explicitly within the puncture-data family 
is also shown for comparison. 
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FIG. 26: The dependence of angular momentum on orbital 
frequency for several approaches. The results indicate that 
the QC sequence and ISCO is entirely consistent with the PN 
treatment at the level of the precision presently possible. 



B. Correspondence with Post- Newtonian 
parameters 

In Ref. 10 Buonanno and Damour (BD) de- 
scribe how to how to use their conservative 2PN order 
Hamiltonian|Q to estimate the location of the last sta- 
ble orbit of a binary black hole system. In that paper 
they go on to evolve such systems at 2.5PN order produc- 
ing an estimate for the gravitational radiation but here 
we are only considering the conservative system. Going 
beyond the 2PN order treatment, Damour, Jaronowski 
and Schafer (DJS) have resummed the conservative part 
of the 3PN Hamilonian to provide a higher-order evalua- 
tion of ISCO ||l|- Until recently, the description of 3PN 
dynamics has been dependent on an unresolved regular- 
ization ambiguity w^, but a recent dimensional regular- 
ization treatment seems to fix this ambiguity with 
the result, ujs = 0- Accepting this result we can ap- 
ply the methods in pl[ to produce a sequence of quasi- 
circular orbits comparable to our QC sequence. There 
are two approaches to PN resummations recommended 
in Ref. |Q . Leaving the reader to find the details in the 
above references we note only that one approach, called 
the — method is modeled after the resummation ap- 
proach of Damour-Iyer-Sathyprakash |43 and the other 
is an extension similar to the Buonanno-Damour treat- 
ment described above. The cleanest way to compare the 
results of these analysis is by looking at the dependence 
of the angular momentum J on the orbital frequency w, 
a gauge invariant comparison. In this view the minimal 
value in the J curves indicates the ISCO. Following DJS 
we have also provided the results of these treatments at 
the IPN and 2PN orders 



Fig. provides a summary of several results. In this 



plot we can assess several important considerations. The 
reliability of the PN treatment is accessible in two ways. 
The degree of convergence with higher PN order provides 
some indication of how much these results can be trusted. 
Earlier hopes that the one-body resummation approach 
might provide accurate waveform information up to and 
even beyond ISCO, possibly minimizing the need for dif- 
ficult numerical treatments are not supported. The one- 
body approach is not converging, in the sense that the 
2PN curve is everywhere closer to the IPN curve than to 
the 3PN results. The f - method seems to show slow 
convergence for well-separated cases with small values of 
u! say 0.05 or less. Likewise, there is notable sensitivity to 
the choice of resummation method for uj < 0.06 or so in 
the 2PN and 3PN cases. For the IPN case, the one-body 
treatment is much closer to the higher order treatments 
than the IPN results. It seems that the careful crafting of 
the one-body resummation method has managed to min- 
imize the significance of 2PN terms for close agreement 
between the IPN and 2PN results, but this treatment 
has not been done as much for the heretofore unknown 
dependence on the higher PN order terms, thus produc- 
ing the appearance of anti-convergence in the one-body 
resummed PN sequence. 

Comparing the post-Newtonian sequences with our QC 
sequence it is notable that the PN results uniformly 
approach the QC sequence with increasing PN order. 
Surely the limit of the PN sequence curves approaches 
something other than the QC curve, but no difference 
between our curve and the apparent limit of the PN se- 
quences is yet discernible at 3PN order. Consistently the 
location of ISCO estimated by the PN treatments ap- 
proaches the parameters of the effective potential ISCO 
data we used to the level of precision at 3PN. This analy- 
sis suggests that the parameters of our ISCO are at least 
as dependable as those given at 3PN and in fact entirely 
consistent with the PN sequence to the degree of pre- 
cision currently possible. The PN results suggest that 
^isco > 0.075 but do not provide a clear upper bound. 

Going beyond the comparison of the invariants J and 
UJ, we can make comparisons in gauge. The ADM gauge 
applied in the one-body approach is a directly general- 
ization of the isotropic coordinates applicable to confor- 
mally fiat data. We can thus compare the separation de- 
pendence directly for the one-body treatments and QC 
sequence. (The gauge transformation for the j^-method 
is impractical.) The comparison in Fig. ^ shows pro- 
found agreement among all treatments all orders for sep- 
arations rADM > 4:Madm- For closer separation the PN 
gauge transformation begins to fail (cannot be solved) in 
some cases. The other relevant quantity is the effective 
potential E = {1 — Mrest/MADM)/i^ shown in Fig. |2^. 
The results are consistent with the conclusions from the 
J-Lu comparison. 

It is worth remarking here that the above agreement is 
also consistent with some recent results using the stan- 
dard 3PN expansion in the standard Taylor form, with- 
out using resummation tecniques p4|. 
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J FIG. 29: Radiated energy for the QC-sequence. 



FIG. 27: The dependence of orbital frequency on ADM-gauge 
separation shows remarkable agreement among all cases. 
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FIG. 28: The dependence of the effective potential Esff on 
orbital frequency for several approaches is consistent with the 
results for angular momentum shown in Fig. Eq 



V. ANALYZING AND TESTING THE FULL 
MODEL 

A fundamental feature of the Lazarus approach for bi- 
nary black hole spacetimes is that have simultaneously 
at hand alternative approaches applicable to the inter- 
face regions of our model spacetime. We have exten- 
sively demonstrated how we can take advantage of this 
feature at the FN-CL interface to cross-check the results 
of the CL and FN treatments against each other, and in- 
deed to better understand the physics of the spacetime. 
Even with our very simple model for the early-part of the 
spacetime we can begin applying the same type of test 
at the FL-FN interface. 



The principle is the same, now, instead of altering the 
time of transition from FN to CL we alter the time of 
transition from FL to FN, and compare the resulting 
waveforms. In practice this means we select alterna- 
tive data sets from the QC sequence and compare the 
waveforms. If our model is astrophysically reliable, then 
the shape of our waveforms should be independent of 
variations in the choice of QC-sequence data set (labeled 
by a QC-number). In this section we will compare the 
waveforms generated over a broad range of the QC se- 
quence, treating the cases ranging from QC-0 to QC-4 
with initial proper separations ranging from about 5.0M 
to about 7.8M. In each case we find a total radiation en- 
ergy of about 2.5-3% M du ring the course of the burst. 
As described in Sec. [II B , we look for a plateau in the 
dependence of the total radiated energy on the FN-CL 

The several 
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transition time T, shown here in Fig. 
methods described in Sec. IIIB collectively suggest lin- 
earization time of around T ~ 9 — 10, 10 — 11, 10 — 11, 11 — 
12, 13 - UM respectively for case QC-0 through QC-4. 
In all cases numerical inaccuracies begin to affect our re- 
sults in the vicinity of T = 12 — 13M as is apparent in 
the figure so the later case, especially QC-4 are only very 
marginally linearized before numerical problems are sig- 
nificant. As is discussed in more detail below, we find 
that significant amounts of angular momentum are radi- 
ated away. For the lower separation cases, QC-0 - QC-2 
roughly .10 — AIM^ of the angular momentum is lost, 
enough to significantly alter the value of the close-limit 
background black hole spin parameter a, so for the re- 
sults presented in this section we have set the background 
black hole angular momentum to the value given in Table 
I minus O.llM^. |^ Likewise we reduce the background 
mass to M = 0.974. 

As our simple early-stage model has less structure than 
the late-stage part, we have a slightly weaker comparison 
to make between the two treatments. 



IIIC 



we 



In Sec. 

had no freedom in our waveform superposition compar- 
ison. Now, for the earlier transition, the shape of the 
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TABLE II: off set shift t ^ {t - to) 
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FIG. 30: Real-part of QC-sequence waveforms for the 1 = 2 
multipole. This plot compares ten waveforms including in 
the FN-CL transition time T as well as variations in the QC- 
sequence from QC-0 to QC-4 (L = 5.0A'/ to L = 7.8Af). The 
agreement over a broad range of model parameters provides 
support for the fidelity of our astrophysical model. 



waveform is still fully constrained, but the overall phase- 
or time-lag between different waveforms is free because 
our simple early-stage model doesn't include any infor- 
mation about the time and rotational offsets between the 
various QC sequence slices as they are rigidly embedded 
in the "astrophysical" spacetime. We will have to set this 
by hand from the waveforms. 



As motivated in Sec. [II C 



we will filter all the wave- 
forms considered here to take only the I = 2 part. This 
has the effect of efficiently eliminating of the higher an- 
gular frequency noise generated in the numerical simula- 
tions, while preserving practically all the physical radi- 
ation. Also, restricting to the I — 2,m — 2 part of the 
radiation as the overwhelmingly predominant component 
completely specifies the angular dependence so that we 
can, without loss of generality, focus on the radiation 
projected toward an observer on the positive z-axis. 

On the basis of the waveform's phase, a good fit results 
from shifting the time-axis by 0.0, -1.6, -2.9, -5.1, and -8.0 
for the QC-0 through QC-4 waveforms, respectively. The 
results are shown in Fig. ^ Table ^ reports those values 
for approximate linearization times corresponding to the 
FN-CL transition. One can interpret the phase-shifts as 
a measure of the differential orbital-plunging times for 
the elements of the QC sequence. 

The remarkable agreement among 'QC waveforms is 
the result of certainly an appropriate combination of or- 
bital parameters X and P. Note that if we keep X fixed 
and change P, like in the 'P-sequence', the waveforms 
do not superpose, and if wc choose similar P and dif- 
ferent initial separations X, like comparing the wave- 
forms of P = 2Pisco/3 and QC-3 (see Figs. |l§]3^) or 
P = 5Pisco/(^ and QC-1 (see Figs. p^[]30| ) we do not find 
any superposition of amplitudes nor phases. 



Name 



L/M 



to/M 



T/M 



QC-0 
QC-1 
QC-2 
QC-3 
QC-4 



5.0 
5.5 
5.9 
6.7 
7.8 



0.0 
1.6 
2.9 
5.1 
8.0 



9- 10 

10- 11 

10- 11 

11- 12 
13 



It is also notable that there seems to be a pretty 
smooth transition from inspiral to plunge. This suggests 
that it will be hard to make a sharp definition of the 
ISCO, but our results compared with those of the Post- 
Newtonian approximation and that of Grand- 
clement et al. ]27|] seems to indicate that starting from 
QC-3 or QC-4 data produces a good approximation to 
the plunge. 

The striking agreement seen in Fig. ^ clearly supports 
the notion that these different initial data sets corre- 
spond approximately to the slices of a single astrophysical 
spacetime. In fact the agreement far exceeds what any- 
one would have expected from this family of conformally 
flat initial data sets, widely regarded on the basis of a pri- 
ori kinematical analysis as needing to be improved upon 
before any astrophysically meaningful results might be 
obtained. Whatever the difficulties are with these data, 
they do not strongly affect the resulting waveforms. 



A. Polarization 

We observe an interesting correlation of the two po- 
larizations of the waveforms here represented by the real 
and imaginary parts of rip 4^ along the 2— axis. Except at 
early times, they are pretty much 90 degrees out of phase. 
The condition corresponds to circular polarization, as 
can be seen more directly if we consider our complex- 
valued waveforms in the representation ^p4{t) = A(t)e*'^*. 
Physically this describes the wave according to an in- 
stantaneous magnitude A and a polarization angle if. 
Fig. ^ shows the time dependence of the polarization 
angle for one of our waveform cases, QC-3. After some 
complicated initial undulations (after around 40M, recall 
that the observer is located at r* « 30M) the polariza- 
tion angle grows smoothly indicating circular polariza- 
tion. Fig. y3 shows the instantaneous magnitude of the 
wave, independent of the polarization angle. For circular 
polarization the magnitude does not vary at the radiation 
frequency, but more slowly. The figure compares wave 
magnitude for the QC-sequence cases. At early times, the 
wave is not circularly polarized, and the magnitude os- 
cillates strongly before becoming dominated by the non- 
oscillating (circularly polarized) component. The curves 
in the figure have been translated by the same time- 
shift applied in Fig. ^ above to get phase agreement in 
the waveforms. The arched circularly-polarized late-time 
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to circulate around the black hole rather than the bell- 
like ringing, one might naively expect. 

We remind the reader that this discussion of waveforms 
is focused on the appearance of the radiation to distant 
observers on the system's rotational axis. For observers 
in other locations the polarization will generally be el- 
liptical, but conforming to a simple pattern consistent 
with the circular motion of the source. In the equato- 
rial plane the observer "sees" no circulation in the source 
and the radiation reduces to polarization. Moving 
around toward the negative z-axis, the observed radiation 
approaches circular polarization of the opposite orienta- 
tion. 



FIG. 31: Polarization angle. For all QC-sequence cases the 
radiation is predominantly circularly polarized after an early 
period of initial-data dependent fluctuations. The colors are 
as labeled in Fig. E9l The curves have been offset if time as 



□ 



B. Radiated Energy and Angular Momentum 

The total radiated energy E and angular momentum 
J at infinity per unit time (u = t — r*) are computed as 
inRef. |8|, 




FIG. 32: Instantaneous radiation magnitude, independent of 
polarization angle. The colors are as labeled in Fig. The 
inset shows the early part with time set to zero at the start 
of numerical evolution. 



parts of these roughly superpose after this time-shift, but 
the oscillatory part, shown without the time offset in the 
inset, is timed to the initial data, and varies notably as 
we change among the different data sets. This suggests 
that the non-circularly-polarized part of the waveform is 
predominantly an initial-data artifact. Notably also, this 
non-astrophysical artifact component seems to shrink for 
the cases of initially more separated systems, as expected 
since the astrophysical interpretation of the initial data is 
not problematic in the limit of well-separated black holes, 
reproducing two Schwarzschild black holes with an small 
boost. 

The physical interpretation of this polarization effect 
is that the system's dynamics are dominated throughout 
the interaction by rotational motion. Even during the 
"ringing" of the final black hole, the perturbations seem 
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(5.2) 



where dQ — sin 9 d9 dip. 

If we consider the magnitude and polarization angle 
representation of h as in Eq. ( |5.3| ), i.e. 2 J J ipi — 
v4(i)e"^'^*-' then we find that radiation energy is gener- 
ated both by time dependence in the magnitude A and 
the polarization angle, but angular momentum flux is 
generated only by the rotation of the polarization vector. 
Thus we can expect to approach maximally efficient ra- 
diation of angular momentum dJ = m/ojdE [p^ for the 
case of circular polarization. This is perhaps the most 
important consequence of the circular polarization pat- 
tern which we observe in our radiation. Note that, for our 
case the most significant radiation has m = 2 azimuthal 
dependence with a frequency near lu ~ 0.5 so that we can 
expect dJ - AMdE. 



C. Memory 

Up to this point we have studied the radiation as repre- 
sented by 1/^4 , which is related to the second time deriva- 
tive of the usual metric disturbance waveform compo- 
nents and /ix by |Q9| 



2 hm 



dt' / dt"'4!4 



(5.3) 




FIG. 33; rh+ and rhx for QC-3 with trivial integration con- 
stants. 

Physically our waveforms correspond to the tidal acceler- 
ation of neighboring observers (or mirrors perhaps). Nat- 
urally we can integrate to get their relative velocities and 
displacements, but this presents us with a problem: How 
do we set the integration constants corresponding to the 
initial displacements and relative velocities of our neigh- 
boring observers? A burst of radiation which allows ini- 
tially quiescent observers to return to their original rest- 
ing places is described as a "normal" burst, otherwise it 
is called a "burst with memory." In trying to integrate 
our waveforms we find that observers at rest when our 
wave arrives will not return to their original resting po- 
sitions, and seem not to return to rest at all. The result 
of direct integration of our waveforms is shown in Fig. ^ 
for a typical case, QC-3. 

The figure shows a strong memory effect for hx and a 
much weaker effect for h+. The hy^ effect is very robust 
with respect to the the QC-sequence, while the mem- 
ory effect may increase slightly with initial separation, 
but this is not clear. Does this mean then that collid- 
ing black holes, contrary to expectation, produce a burst 
with memory? Not necessarily. Our wave is only the 
last part of a much longer wavetrain. Let us adopt the 
hypothesis that the full wave train is a "normal" burst, 
and consider whether our results contradict the hypoth- 
esis. A key point is that neighboring observers of the 
full wavetrain, which were initially quiescent, will not be 
at rest in the middle of the wavetrain, where our wave 
starts. By hypothesis, we expect these observers to re- 
turn to relative rest so we must expect a change in their 
relative velocities and positions through the passing of 
our part of the wave. Our results are consistent with 
the hypothesis of finally resting observers if the "cross" 
observers, separated transversely by AL at a distance R 
from the binary system, are initially moving with speed 
V ^ 0.0014AL /R and the "plus" observers are initially 
displaced by about 0.02 AL M/ R. From the quadrupole 
formula for two orbiting point masses positioned instan- 
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FIG. 34: r/i-|_ and rhx with integration constants set to fit 
"normal burst hypothesis". Solid curves are '4-' polarization. 
Dotted curves are cross polarization. 



taneously on the y-axis as our black holes are initially we 
expect that the cross observers will be moving with an 
initial velocity of 2QA and the plus observers will be ini- 
tially displaced by A where is the orbital velocity and 
A is the momentary amplitude of the radiation. Our re- 
sults are then consistent, at an order of magnitude level 
with what might be expected from extrapolated post- 
Newtonian results as presented in At present it 
would be inappropriate to pursue this quantitative rela- 
tionship too far since there may be subtle effects, such as 
transients in the early parts of our waveforms or frame 
dragging of the radiation. In the end, our memory effect 
is plausibly consistent with the "normal" burst hypoth- 
esis. But this is certainly an interesting area for further 
study. 

If we provisionally adopt the normal burst hypothesis, 
it is then most reasonable to perform the integrations 
required in Figure ^ with the integration constants set 
such that the strains approach zero at the end of the 
burst. Because our waves are not completely quiescent 
and prone to small errors in the late time region, some 
visual judgment is required in choosing appropriate inte- 
gration constants. After shifting the time as above, we 
get reasonable waveforms upon offsetting the /i+ curves 
by 0.015, 0.017, 0.020 and 0.027 respectively for the cases 
QC-0 through QC-3 in the second integral only. For hx 
we offset by - 0.002 x t - 0.009 The resulting 

waveforms are shown in Fig. Unlike the other plots 
in this section we have included only one curve for each 
polarization from each QC data set. It is interesting to 
compare these "strain" waveforms with the other "accel- 
eration" waveforms presented so far in this paper. An 
obvious difference is that the sharp growth in amplitude 
during the burst, seen in Fig. |3^ for instance, is greatly 
muted here. The stronger growth in the "04 waveforms is 
a consequence of the two factors of, uj, which is rapidly 
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TABLE III: 3PN-E1B-ISCO data 



Name 


L/M 


±X/M 


±P/M 




Mlu m/M 


3PN-E1B 


~ 7.3 


2.119 


0.200 


0.8476 


~ 0.085 0.469 



increasing, obtained from two time derivatives. A con- 
sequence of this evident in the figure is that the peak in 
amphtude comes earher in the strain waveforms, pushing 
back to around, T = 40M near the beginning of the as- 
trophysically credible part of the waveform. In the early 
time region the figure suggests that the large QC cases 
are tending toward longer wavelength early waveform, as 
one would hope to see when beginning with "astrophys- 
ically earlier" data. Finally we note that phase relation 
between the two polarizations, which we have attributed 
to circular polarization is not disrupted by two time inte- 
grations, and is equally evident in the strain waveforms. 
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FIG. 35: Comparisons of Waveforms 3PN, QC-3, QC-4 Wave- 
forms 



D. Comparison with Post- Newtonian results 



VI. SUMMARY OF ASTROPHYSICAL 
RESULTS 



As a point of comparison we also consider the 3PN 
ISCO determined by the method of the effective one 
body (ElB) of Ref. |42j and obtain the parameters (in 
the ADM gauge) given in table [II. Note that in Ref. 
p4 1 it was stressed the fact that to third post- Newtonian 
order is not obvious that the one body resummation is 
better than the 'bare' results. In this case we find that 
the ISCO so determined have parameter very close to the 
QC — 2 case studied here. 

In order to quantify the gravitational radiation gener- 
ated during this plunge we use this 3PN ISCO parameters 
to evolve Bowen-York initial data. The resulting wave- 
form is plotted in Fig. ^ and compared with the QC-3 
and QC-4 cases, showing closer agreement with the QC-3 
waveform. 

For the full nonlinear part of the simulation we have 
used a nonuniform grid of 512^ x 256 and 448^ x 224 
with boundaries at 37M We have also considered differ- 
ent nonuniform grids and boundary locations on another 
512^ X 256 simulation to check the robustness of the re- 
sults. The largest runs lasted for 10 clock wall hours on 
the 64 nodes of the SR-8000 Hitachi computer at the LRZ 
and used up to 170Gb of RAM memory. 

The analysis of the FN determination of the ISCO 
shows that its separation and linear momentum parame- 
ters are close to those of the case QC-3. In retrospective, 
as we have seen, the plunge waveform does not dramati- 
cally change for the five cases studies, i.e. QC-0 to QC- 
4, hence our predictions for the plunge waveforms are 
compatible with the third Post-Newtonian order orbital 
parameters for the ISCO as well. 



In the last section we have presented a detailed anal- 
ysis of resulting waveforms from our model for the fi- 
nal plunge of a system of equal-mass, non-spinning (ir- 
rotational), non-elliptically orbiting binary black holes, 
results. The domain of astrophysical validity waveform 
is limited to the time from right around the first burst 
peak forward but the surprisingly lack of internal model- 
dependence in our results suggests that we can be con- 
fident (say to the level of 20%) in our description of the 
final radiation burst from such a system. We find these re- 
sults convincing enough to merit their adoption, among 
interested researchers, as a provisional description for rel- 
evant astrophysical work, such as in the development of 
strategies for observation of gravitational waves. As it 
has been emphasized that any new information may be 
very valuable for improving the cfhciency of gravitational 
wave search algorithms [Q , and for the interest of readers 
not yet immersed in the effort to model relativistic sys- 
tems we briefly review the results from a more applied 
astrophysical viewpoint. 

Our waveforms cover the final few cycles of gravita- 
tional radiation from the peak onward. Within about 
four cycles the signal is reduced to only a few percent 
of its peak value. In geometric units this takes about 
50Af, which corresponds, for a 2OM0 system to about 
0.1ms or, for a lO^M© system of supermassive black holes 
about one minute. During this brief time vast amounts 
of energy are released, equivalent to about 2.5-3% of the 
system's total mass and ~ 12% of its angular momen- 
tum. The peak gravitational wave luminosity is about 
0.0015c^/G or about 5 x lO^^ergs/s. The radiation is 
strongly polarized in the same manner as expected from 



a rotating gravitational source with circular polarization 
for radiation propagating along the rotational axes and 
"+" polarization for radiation in the equatorial plane. As 
generally expected the radiation is predominated by an 
I = 2, m = 2, s = —2 spin- weighted spherical harmonic 
angular distribution so that that the intensity is 16 times 
greater along the axes than in the equatorial plane. 

The simple appearance of our waveforms, as shown in 
Fig. ^ suggests that rather than asking for awkward nu- 
merical data files, other researchers may benefit from a 
straightforward analytic fitting of our result with essen- 
tially the full information at our level of confidence. For 
a representative case, the QC-3 waveform, we can impose 
circular polarization, a frequency and amplitude depen- 
dence of the form 
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A simple interpretation of the parametrization we 
have chosen for the plunge waveform starts from the 
amplitude-phase representation of the waveform r'ip4 as 
above. We then model the time dependence of the orbital 
frequency of the last stage of the inspiral by a linear in- 
crease that changes slope through the plunge piece of 
the waveform and then matched the final quasi-normal 
frequency of the ringing down remnant black hole. For 
the time variation of the amplitude we have chosen to 
match between an initial exponential growth and a final 
exponential decay. 

The quantities labeled "QNM" are set to their ex- 
pected values for the least damped quasi-normal mode 
from black hole perturbation theory, in this case wqnm — 
0.55 and oqnm = —0.073. The quantity ao is chosen in- 
formally to have the value 0.0085. The other eight quan- 
tities are chosen by non-linear least-squares fit to have 
the approximate respective values {ujQ — 0.2894, t^o — 
33.00, t^i = 67.05, ctO = 0.2192, t^o = 32.54, t^i = 
62.92, (po — —4.40, oq — —6.31}. A comparison of the 
fit to the original curve is shown in Fig. |3^. We plan 
to further investigate the "fitting" representation of the 
waveforms to optimize the number and robustness of the 
parameters to be used as well as the choice of the fitting 
functions. This will also facilitate the matching of the 
plunge waveforms with a, for instance, Post-Newtonian 
one for the inspiral phase. 

The curve shown for rip/^{t) displays only the time de- 



1 ' 1 ' 


1 1 , 1 , 1 , 






.... f„ 








.... fit 




















III, 


1 < 1 < 1 < 1 



m 



FIG. 36: Comparison with the fitted curve. 



pendence. For the full wave information, we have 



Mt-r\e,^) = -2Y22^t~^^\t-r*) 

+-2Y2-2 4"'=-^\t-r*) (6. 

where the sYim are spin-weighted spherical harmonics 



-2Y22 = V574^cos4(^/2)exp(2i(^) (6.7) 
V5/47r sin'' (61/2) exp(-2i(p). (6.8) 
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VII. DISCUSSION 

We have produced the first astrophysical model for 
binary black hole coalescence waveforms. We use the 
Lazarus method applying the techniques of numerical rel- 
ativity in the strongly interacting intermediate plunge re- 
gion and close-limit black hole perturbation theory in the 
late-time region. As a stand-in for the "far-limit" early 
part of the dynamics we have taken data from a family 
of solutions to Einstein's equations restricted to a space- 
like slice of spacetime. This Lazarus approach is ideal in 
this application for which each component treatment is 
at best marginally applicable by itself. Not only does the 
Lazarus method expand the scope of treatable problems, 
but it provides a natural framework for cross-checking 
the various treatments against each other. 

In the first part of the paper we have performed a thor- 
ough study of the performance the late-time treatment, 
including the numerical relativity and close-limit regions 
on a series of initial data set approaching the fiducial 
target "ISCO" data, the results of which were recently 
reported on in Ref. the present study strongly es- 

tablishes the effectiveness of our late-time methods. 

Our early-time model is based on a selection from one 
of the standard families of black hole initial data used in 
numerical relativity, the QC-sequence. Because of the 
generic uncertainty in describing astrophysically accu- 
rate initial data consistent with Einstein's equations for 
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strongly interacting black holes, we had little basis, a pri- 
ori, for trusting these data as astrophysically reasonable. 
Thus, the main part of this paper has addressed the ro- 
bustness of the overall results against variations in the 
physical transition point between the early and late-time 
portions of our model. In other words, we have tested the 
kinematically defined initial data sets against numerical 
relativistic dynamics. The results show an overwhelming 
correspondence, which strongly suggests that, at least 
for, non-spinning black holes, the subtleties of selecting 
initial data are not as consequential to the resulting wave- 
forms as has been generally expected in the community, 
and supporting our results as a rough first-look at the 
gravitational radiation which can be expected from the 
coalescence of equal-mass non-spinning black holes. 

We find the radiation to be dominated by circularly 
polarized, I = 2, \m\ = 2 waves radiating 3% of the sys- 
tem's total energy in a few cycles and significant angular 
momentum, so that the remnant black hole has an an- 
gular momentum parameter of roughly a ^ 0.7M, as we 
have summarized the results in relation to astrophysi- 
cal gravitational wave observation observation in the last 
section. 

In the introductory remarks, we noted the importance 
of having at least a provisional astrophysical coalescence 
study not only for the results it provides to observers, 
who need information for developing their observational 
strategies, but also for the benefit of theorists who wish 
to direct their effort efficiently toward the purpose of as- 
trophysical application. In this light, we conclude with a 
characterization, suggested by the results, of the course 
we will pursue in future work. We first note that past 
developments in numerical relativity have been critical 
in our application and have already brought the state 
of the art to the point where it can be applied astro- 
physically. One of the great theoretical worries has been 
that astrophysical initial data for relatively close black 
holes might require a great deal of further work. Our 
results suggest this problem may not be critical, but fur- 
ther work is needed when the individual spins of the in- 



teracting black holes may become relevant. It will also be 
interesting to learn how generally this waveform robust- 
ness extends when studying other initial data families. A 
weakness of our present early-time model is that it pro- 
vides no direct connection to the inspiral radiation, and 
thus no cross check with post-Newtonian (PN) dynamics. 
Further work in this direction will benefit from work on 
initial data which includes PN dynamical radiation infor- 
mation, as well as longer-running numerical simulations 
allowing us to extend waveforms farther back in time. 
Considering the waveforms here presented, we are im- 
pressed by its simplicity, a simple multipolar description 
with monotonically increasing frequency and smoothly 
changing amplitude. However, for systems of strongly 
spinning black holes or with significantly unequal masses 
other arguments indicate complicated waveforms result- 
ing from very involved interactions. A significant ob- 
servationally motivated goal will then be to learn how 
generally, in the parameter space of plausible black hole 
pairs, we can expect the qualitatively simple type inter- 
action and radiation. We expect such a characterization 
to the parameter space to require significant additional 
effort, but we can begin soon by comparatively studying 
co-rotational systems and exploring the lowest order spin 
and unequal mass effects more quickly. 
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